Statistique exploratoire spatiale : Travaux Pratiques N°11

ENSAE Pierre NDIAYE de Dakar ISE1-Cycle long 2024-2025

COURS DE Statistique exploratoire spaciale avec M.Aboubacar HEMA

Travaux Pratiques N°11

Groupe : Logiciel R

Pays : Mali

Membres : Fogwoung Djoufack Sarah-Laure, SENE Malick, Niass Ahmadou, Celina Nguemfouo

===================== CONSIGNE =================

1. Faire au préalable tout ce qui a été fait au TP6
2. Binariser le raster des évènements de telle sorte que les pixels ayant plus de 5 evenements prennent la valeur 1 et dans le cas contraire, c’est la valeur 0 qui est affectée
3. Pour le raster population, amener la résolution qui était de 100m à 5km par agrégation
4. Calculer le confliction Diffusion Indicator, au niveau pays et au niveau région, en procédent comme suit: a- Binariser le raster population de 5km, qui prend ainsi la valeur 1 si le nombre d’habitants est supérieur à 50 et 0 sinon b- Multiplier les deux raster binarisés c- Compter le nombre de 1 dans le raster binarisé de 5km de la population d- Faites la meme chose pour le raster produit obtenu à l’étape b e- Faites le rapport de ce qui est obtenu au d par ce qui est obtenu au c et ce rapport est le CDI recherché

STEP 1: Chargement des bibliothèques nécessaires

suppressMessages({
# Charger les bibliothèques nécessaires
library(ggplot2)  # Pour créer des visualisations statiques.
library(dplyr)    # Manipuler et transformer des données d'un tableau
library(leaflet)  # Créer des cartes interactives.
library(sf)       # Pour les données spatiales
library(raster)   # Manipuler et analyser des données raster
library(terra)    # Pour la gestion des données raster et vectorielles, similaire à raster, mais plus rapide pour les grandes données et avec des fonctionnalités supplémentaires
library(leaflet.extras) # Ajouter des fonctionnalités supplémentaires à leaflet (fournir des controles, recentrage )
library(viridis)  # Générer des palettes de couleurs
library(exactextractr)  #   Calculer des statistiques sur des raster 
library(ggspatial)
})

## Warning: le package 'ggplot2' a été compilé avec la version R 4.4.2

## Warning: le package 'dplyr' a été compilé avec la version R 4.4.2

## Warning: le package 'leaflet' a été compilé avec la version R 4.4.2

## Warning: le package 'sf' a été compilé avec la version R 4.4.2

## Warning: le package 'raster' a été compilé avec la version R 4.4.2

## Warning: le package 'sp' a été compilé avec la version R 4.4.2

## Warning: le package 'terra' a été compilé avec la version R 4.4.2

## Warning: le package 'leaflet.extras' a été compilé avec la version R 4.4.2

## Warning: le package 'viridis' a été compilé avec la version R 4.4.2

## Warning: le package 'viridisLite' a été compilé avec la version R 4.4.2

## Warning: le package 'exactextractr' a été compilé avec la version R 4.4.2

## Warning: le package 'ggspatial' a été compilé avec la version R 4.4.2

Definition du repertoire de travail

setwd("C:/Users/DELL/Desktop/Célina❤/ISEP/ISEP3 2024-2025/Semestre 1/Statistique exploratoire spatiale/TP11")

REVISION DU TP6

Importation des données depuis le fichier CSV envoyé

data <- read.csv("Points_data.csv")

Chargement du fichier de données spatiales (shapefile) pour les limites administratives du Mali

shp <- st_read("DONNEES_MALI/mli_admbnda_adm1_1m_gov_20211220.shp", quiet= TRUE)

Quelques informations sur la base de données (type de l’objet, noms des colonnes, types de données, premieres valeurs)

str(data)

## 'data.frame':    87225 obs. of  31 variables:
##  $ event_id_cnty     : chr  "MLI33140" "BFO12659" "BFO12661" "BFO12665" ...
##  $ event_date        : chr  "04 October 2024" "04 October 2024" "04 October 2024" "04 October 2024" ...
##  $ year              : int  2024 2024 2024 2024 2024 2024 2024 2024 2024 2024 ...
##  $ time_precision    : int  1 1 1 1 1 1 1 1 2 1 ...
##  $ disorder_type     : chr  "Political violence" "Political violence" "Political violence" "Political violence" ...
##  $ event_type        : chr  "Explosions/Remote violence" "Explosions/Remote violence" "Battles" "Battles" ...
##  $ sub_event_type    : chr  "Air/drone strike" "Air/drone strike" "Armed clash" "Armed clash" ...
##  $ actor1            : chr  "Military Forces of Mali (2021-)" "Military Forces of Burkina Faso (2022-)" "JNIM: Group for Support of Islam and Muslims" "JNIM: Group for Support of Islam and Muslims" ...
##  $ assoc_actor_1     : chr  "" "" "" "" ...
##  $ inter1            : chr  "State forces" "State forces" "Rebel group" "Rebel group" ...
##  $ actor2            : chr  "Civilians (Niger)" "JNIM: Group for Support of Islam and Muslims" "Military Forces of Burkina Faso (2022-)" "VDP: Volunteer for Defense of Homeland" ...
##  $ assoc_actor_2     : chr  "Miners (Niger)" "" "VDP: Volunteer for Defense of Homeland" "" ...
##  $ inter2            : chr  "Civilians" "Rebel group" "State forces" "Identity militia" ...
##  $ interaction       : chr  "State forces-Civilians" "State forces-Rebel group" "State forces-Rebel group" "Rebel group-Identity militia" ...
##  $ civilian_targeting: chr  "Civilian targeting" "" "" "" ...
##  $ iso               : int  466 854 854 854 854 854 288 466 566 566 ...
##  $ region            : chr  "Western Africa" "Western Africa" "Western Africa" "Western Africa" ...
##  $ country           : chr  "Mali" "Burkina Faso" "Burkina Faso" "Burkina Faso" ...
##  $ admin1            : chr  "Kidal" "Cascades" "Sahel" "Nord" ...
##  $ admin2            : chr  "Abeibara" "Comoe" "Soum" "Yatenga" ...
##  $ admin3            : chr  "Tinzawatene" "Ouo" "Djibo" "Ouahigouya" ...
##  $ location          : chr  "Zakak" "Dida Forest" "Djibo" "Aorema" ...
##  $ latitude          : num  19.6 10 14.1 13.7 12.7 ...
##  $ longitude         : num  2.891 -4.029 -1.642 -2.334 -0.131 ...
##  $ geo_precision     : int  1 2 2 1 2 1 1 1 2 1 ...
##  $ source            : chr  "FAMAMali; Twitter; Undisclosed Source" "Undisclosed Source" "Al Zallaqa; Whatsapp" "Whatsapp" ...
##  $ source_scale      : chr  "Local partner-Other" "Local partner-Other" "New media" "New media" ...
##  $ notes             : chr  "On 4 October 2024, the Malian air force carried out an airstrike against a convoy in the village of Zakak (Abei"| __truncated__ "On 4 October 2024, the Burkinabe force carried airstrikes against JNIM militants in the Dida Forest (Ouo, Comoe"| __truncated__ "On 4 October 2024, JNIM militants carried out several attacks against positions of Burkinabe force and voluntee"| __truncated__ "On 4 October 2024, JNIM militants attacked a position of volunteer fighters (VDP) in the village of Aorema (Oua"| __truncated__ ...
##  $ fatalities        : int  7 5 0 0 0 0 0 3 0 2 ...
##  $ tags              : chr  "" "" "" "" ...
##  $ timestamp         : int  1728335020 1728358478 1728358478 1728358478 1728358478 1728358478 1728358486 1728358502 1728358504 1728358504 ...

Visualisation des noms des variables

colnames(data)

##  [1] "event_id_cnty"      "event_date"         "year"              
##  [4] "time_precision"     "disorder_type"      "event_type"        
##  [7] "sub_event_type"     "actor1"             "assoc_actor_1"     
## [10] "inter1"             "actor2"             "assoc_actor_2"     
## [13] "inter2"             "interaction"        "civilian_targeting"
## [16] "iso"                "region"             "country"           
## [19] "admin1"             "admin2"             "admin3"            
## [22] "location"           "latitude"           "longitude"         
## [25] "geo_precision"      "source"             "source_scale"      
## [28] "notes"              "fatalities"         "tags"              
## [31] "timestamp"

Differents pays dans la base

unique(data$country)

##  [1] "Mali"          "Burkina Faso"  "Ghana"         "Nigeria"      
##  [5] "Benin"         "Guinea"        "Senegal"       "Ivory Coast"  
##  [9] "Guinea-Bissau" "Mauritania"    "Niger"         "Togo"         
## [13] "Liberia"       "Cape Verde"    "Gambia"        "Sierra Leone"

Conversion des données en un objet spatial sf

data_spatial <- sf::st_as_sf(data, coords = c("longitude", "latitude"), crs = st_crs(shp))

• coords spécifie les colonnes contenant les coordonnées (longitude et latitude)
• crs définit le système de coordonnées à associer aux points, basé sur le shapefile

Visualisation des données spatiales

Création d’une carte statique avec ggplot2

ggplot(data_spatial) +
  geom_sf(fill = NA, color = "blue", size = 0.5) +
  aes(colour = country) + # Pour colorer les points selon le pays
  geom_sf(size = 1.2) +   # Ajuster la taille des points
  scale_fill_hue(direction = 1) + # Palette de couleurs pour le remplissage
  scale_color_hue(direction = 1)  # Palette de couleurs pour les contours

Carte interactive

## Palette de couleurs
country_palette <- colorFactor(
  viridis(length(unique(data$country)),option = "turbo"), 
  domain = data$country)

## Générer une carte interactive 
leaflet() %>%
  addTiles() %>%  # Ajouter une couche de base (OpenStreetMap)
  # Ajouter les limites administratives depuis le shapefile
  addPolygons(data = shp, color = "brown", weight = 2, opacity = 1, fillOpacity = 0.5,
              popup = ~ADM1_FR) %>%  # Afficher l'information dans une popup
  
  # Ajouter les points d'événements géolocalisés
  addCircleMarkers(data = data_spatial, weight = 0.1, opacity = 2, fillOpacity = 1.4,
                   radius = 1.5,  # Taille des cercles
                   color = ~country_palette(country)) %>% # Couleur selon le pays
  addLegend(    # Ajouter une légende pour indiquer la correspondance des couleurs avec les pays
    "bottomright", pal = country_palette, values = data$country,
    title = "Evenements par pays", opacity = 1
    ) %>%
  addResetMapButton()%>%  # Recentrer la carte
  addFullscreenControl()  # Contrôle pour passer en mode plein écran

Sélection du Mali

              #### Specifying our Area of interest (AOI)
    # On selectionne donc le pays
AOI = "Mali"

    # Créer un sous-ensemble de données ne contenant que les événements relatifs au Mali
AOI_event <- data_spatial %>%
  filter(country == AOI)

    # Visualisation des événements au Mali
ggplot(AOI_event) +
  aes(fill = event_type, colour = event_type) + # Associer les couleurs aux types d'événements
  geom_sf(size = 1.2) + # Représentation des géométries spatiales
  scale_fill_hue(direction = 1) + # Échelle de couleur pour les types d'événements
  theme_minimal()

## Palette de couleurs
event_palette <- colorFactor(palette = "Set2", domain = AOI_event$event_type)

    # Créer une carte interactive
leaflet() %>%
  addTiles() %>%  # Couche OpenStreetMap
  # Ajouter les limites  (administration de niveau 0 - national)
  addPolygons(data = shp, color = "brown", weight = 2, opacity = 1, fillOpacity = 0.5,
              popup = ~ADM1_FR) %>%  # Afficher l'information dans une popup
  
  # Ajouter les points d'événements 
  addCircleMarkers(data = AOI_event, weight = 0.1, opacity = 2, fillOpacity = 1.4,
                   radius = 2,  # Adjust circle size
                   color = ~event_palette(event_type)) %>%
  addLegend("bottomright", pal = event_palette, values = AOI_event$event_type,
            title = "Event Type", opacity = 1) %>%
  addResetMapButton()%>%  # Recentrer la carte
  addFullscreenControl()  #ajout du basculement en mode plein écran

Calcul du nombre d’évènements par admin (0-3)

Association des points d’événements (data_spatial) aux limites administratives données par le shapefile(shp)

points_ml<- st_join(data_spatial, shp, join = st_intersects)

points_ml %>% data.frame() %>% tail(5)# Afficher les 5 dernières lignes du DataFrame résultant pour vérifier la jointure

##       event_id_cnty      event_date year time_precision          disorder_type
## 87221          SIE2 01 January 1997 1997              3     Political violence
## 87222          SIE3 01 January 1997 1997              3     Political violence
## 87223          SIE6 01 January 1997 1997              3     Political violence
## 87224          SIE7 01 January 1997 1997              3     Political violence
## 87225          SIE8 01 January 1997 1997              3 Strategic developments
##                       event_type                   sub_event_type
## 87221                    Battles     Government regains territory
## 87222                    Battles                      Armed clash
## 87223 Violence against civilians                           Attack
## 87224                    Battles                      Armed clash
## 87225     Strategic developments Headquarters or base established
##                                            actor1 assoc_actor_1
## 87221 Military Forces of Sierra Leone (1996-1997)              
## 87222                             Kamajor Militia              
## 87223 Military Forces of Sierra Leone (1996-1997)              
## 87224                             Kamajor Militia              
## 87225 Military Forces of Sierra Leone (1996-1997)              
##                  inter1                                      actor2
## 87221      State forces             RUF: Revolutionary United Front
## 87222 Political militia             RUF: Revolutionary United Front
## 87223      State forces                    Civilians (Sierra Leone)
## 87224 Political militia Military Forces of Sierra Leone (1996-1997)
## 87225      State forces                                            
##       assoc_actor_2       inter2                    interaction
## 87221                Rebel group       State forces-Rebel group
## 87222                Rebel group  Rebel group-Political militia
## 87223                  Civilians         State forces-Civilians
## 87224               State forces State forces-Political militia
## 87225                                         State forces only
##       civilian_targeting iso         region      country   admin1  admin2
## 87221                    694 Western Africa Sierra Leone  Eastern    Kono
## 87222                    694 Western Africa Sierra Leone Southern      Bo
## 87223 Civilian targeting 694 Western Africa Sierra Leone Southern  Bonthe
## 87224                    694 Western Africa Sierra Leone Southern Moyamba
## 87225                    694 Western Africa Sierra Leone Southern Moyamba
##        admin3    location geo_precision                           source
## 87221   Gbane       Mandu             2 No Peace Without Justice; SL-LED
## 87222 Selenga     Selenga             2 No Peace Without Justice; SL-LED
## 87223  Imperi York Island             2 No Peace Without Justice; SL-LED
## 87224 Fakunya     Fakunya             2 SL-LED; No Peace Without Justice
## 87225    Kori      Taiama             2 No Peace Without Justice; SL-LED
##                  source_scale
## 87221 Local partner-New media
## 87222 Local partner-New media
## 87223 Local partner-New media
## 87224 Local partner-New media
## 87225 Local partner-New media
##                                                                                                                                                                               notes
## 87221 Around 1 January 1997 (month of), Military Forces of Sierra Leone (1996-1997) and RUF: Revolutionary United Front engaged in battles in Mandu (Eastern, Kono). No fatalities.
## 87222                            Around 1 January 1997 (month of), Kamajor Militia and RUF: Revolutionary United Front engaged in battles in Selenga (Southern, Bo). No fatalities.
## 87223         Around 1 January 1997 (month of), Military Forces of Sierra Leone (1996-1997) engaged in violence against civilians in York Island (Southern, Bonthe). No fatalities.
## 87224           Around 1 January 1997 (month of), Kamajor Militia and Military Forces of Sierra Leone (1996-1997) engaged in battles in Fakunya (Southern, Moyamba). No fatalities.
## 87225                                                                                                                                                                          base
##       fatalities tags  timestamp Shape_Leng Shape_Area ADM1_FR ADM1_PCODE
## 87221          0      1670286851         NA         NA    <NA>       <NA>
## 87222          0      1670286851         NA         NA    <NA>       <NA>
## 87223          0      1670286851         NA         NA    <NA>       <NA>
## 87224          0      1670286851         NA         NA    <NA>       <NA>
## 87225          0      1678830926         NA         NA    <NA>       <NA>
##       ADM1_REF ADM1ALT1FR ADM1ALT2FR ADM0_FR ADM0_PCODE date validOn validTo
## 87221     <NA>       <NA>       <NA>    <NA>       <NA> <NA>    <NA>    <NA>
## 87222     <NA>       <NA>       <NA>    <NA>       <NA> <NA>    <NA>    <NA>
## 87223     <NA>       <NA>       <NA>    <NA>       <NA> <NA>    <NA>    <NA>
## 87224     <NA>       <NA>       <NA>    <NA>       <NA> <NA>    <NA>    <NA>
## 87225     <NA>       <NA>       <NA>    <NA>       <NA> <NA>    <NA>    <NA>
##                      geometry
## 87221 POINT (-10.9332 8.4642)
## 87222 POINT (-11.7047 8.1221)
## 87223 POINT (-12.4694 7.5317)
## 87224   POINT (-12.338 8.231)
## 87225   POINT (-12.06 8.2013)

Garder les points du Mali

points_ml <- points_ml %>% filter(!is.na(ADM1_PCODE)) # Supprimer les points qui ne sont pas associés à une région administrative (valeurs NA dans ADM1_PCODE)
points_ml %>% data.frame() %>% dim() # 11547

## [1] 11547    42

Donc 6 points à l’extérieur du raster

dim(AOI_event) # 11541

## [1] 11541    30

Nombre de points (événements) pour chaque région administrative (ADM1_FR)

point_counts <- points_ml %>%
  group_by(ADM1_FR) %>%  #  groupper par régions
  summarise(Nombre_attaques = n())
point_counts %>%
  st_drop_geometry() %>%
  data.frame() # Affichage dans un tableau

##       ADM1_FR Nombre_attaques
## 1      Bamako             446
## 2         Gao            1956
## 3       Kayes             265
## 4       Kidal             935
## 5   Koulikoro             429
## 6      Menaka             720
## 7       Mopti            3852
## 8     Sikasso             308
## 9       Ségou            1380
## 10 Tombouctou            1256

Tableau croisé : Nombre d’événements par région administrative (admin1) et type d’événement

t1 <-table(AOI_event$admin1, AOI_event$event_type) 
t1

##             
##              Battles Explosions/Remote violence Protests Riots
##   Bamako          31                          6      245    65
##   Gao            476                        251       78    30
##   Kayes           79                          4       75    28
##   Kidal          237                        368       44    27
##   Koulikoro      132                         54       40    26
##   Menaka         211                         93       13     7
##   Mopti         1095                        521      105    39
##   Segou          410                        199       35    14
##   Sikasso         79                         25       44    44
##   Tombouctou     300                        216       59    44
##             
##              Strategic developments Violence against civilians
##   Bamako                         58                         41
##   Gao                           444                        679
##   Kayes                          28                         51
##   Kidal                         135                        124
##   Koulikoro                      72                        103
##   Menaka                        137                        246
##   Mopti                         599                       1498
##   Segou                         231                        494
##   Sikasso                        58                         58
##   Tombouctou                    225                        411

Nombre d’événements par division administrative (admin2)

t2 <- table(AOI_event$admin2)  %>% data.frame() 
colnames(t2) <- c("division", "nbre_evenements")
head(t2,10)

##          division nbre_evenements
## 1        Abeibara             117
## 2  Anderamboukane             109
## 3         Ansongo             824
## 4       Bafoulabe              32
## 5          Bamako             446
## 6         Banamba              52
## 7      Bandiagara             763
## 8         Bankass             481
## 9       Baraoueli              12
## 10            Bla              14

Fonction pour créer un raster basé sur les données géolocalisées

# Fonction pour créer un raster basé sur les données géolocalisées
Create_raster <- function(datafile, filename ="Rasterisation.tif") {
  # Reprojection des données spatiales pour un système de coordonnées en mètres
# EPSG 32629 : UTM Zone 29N, adapté au Mali
  AOI_prj <- st_transform(datafile, crs = 32629)  
  
# Définir l'étendue géographique (extent) à partir des limites de l'objet spatial reprojeté
  ext <- raster::extent(sf::st_bbox(AOI_prj))  # Conversion des limites en format `extent
# Spécification de la résolution en mètres (5 km ici)
  res <- 5000  
# Définir un système de coordonnées pour le raster 
  rast_crs <- CRS("+proj=utm +zone=29 +datum=WGS84 +units=m +no_defs")
  raster_template <- raster::raster(ext=ext, resolution=res, crs=rast_crs)
# Rasteriser les données : calculer la somme des événements dans chaque cellule de la grille
Raster <- raster::rasterize(AOI_prj,raster_template,field=1, fun= sum)
# Sauvegarder le raster en format GeoTIFF
  raster::writeRaster(Raster, filename = filename, format = "GTiff", overwrite = TRUE)
  return(Raster) # Pour retourner l'objet raster
}

Application de la fonction sur les données des événements

AOI_Raster <- Create_raster(AOI_event, "Rasterisation_general.tif")

Fonction pour afficher une carte

Create_map <- function(raster){
  ## Palette de couleurs
  pal <- colorNumeric(palette = viridis(1000, option = "viridis"), 
                      domain = raster::values(raster),
                      na.color = "transparent")
# Créer une carte interactive
  leaflet() %>%
    addTiles() %>%  # Couche de base (OpenStreetMap)
    
    addPolygons(data = shp, color = "brown", weight = 2, opacity = 0.2, fillOpacity = 0.1,
                popup = ~ADM1_FR) %>%  
    
    addRasterImage(raster, opacity = 2,colors= pal) %>%
    addLegend(pal = pal, values = raster::values(raster),
              title = "Nombre d'événement") %>%
    addResetMapButton() %>%  # Recentrer la carte
    addFullscreenControl()  #ajout du basculement en mode plein écran  
}

Application de la fonction sur le raster

Create_map(AOI_Raster) 

## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Résumé statistique des valeurs du raster

summary(raster::values(AOI_Raster)) 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     1.0     1.0     2.0     5.8     4.5   445.0   94959

Création de raster pour chaque année et visualistion du nombre d’attaques par année

for( i in unique(data$year)){
  assign(paste0("data_", i), data[data$year == i,])
  assign(paste0("AOI_Raster_", i), Create_raster(AOI_event,  paste0("Rasterisation_", i, ".tif")))# Création du raster pour chaque année et sauvegarde en fichier GeoTIFF
}

Deux premieres lignes

head(data_2020, 2)

##       event_id_cnty       event_date year time_precision      disorder_type
## 42787       CDI2482 31 December 2020 2020              1 Political violence
## 42788      NIG19389 31 December 2020 2020              1     Demonstrations
##                       event_type        sub_event_type
## 42787 Violence against civilians                Attack
## 42788                      Riots Violent demonstration
##                                       actor1 assoc_actor_1            inter1
## 42787 Unidentified Armed Group (Ivory Coast)               Political militia
## 42788                      Rioters (Nigeria)                         Rioters
##                        actor2                         assoc_actor_2    inter2
## 42787 Civilians (Ivory Coast) Government of the Ivory Coast (2011-) Civilians
## 42788                                                                        
##                       interaction civilian_targeting iso         region
## 42787 Political militia-Civilians Civilian targeting 384 Western Africa
## 42788                Rioters only                    566 Western Africa
##           country admin1           admin2     admin3   location latitude
## 42787 Ivory Coast  Comoe Indenie-Djuablin Abengourou Abengourou   6.7297
## 42788     Nigeria   Ogun     Egbado South              Oke Odan   6.7000
##       longitude geo_precision                    source           source_scale
## 42787   -3.4964             2 Afrique sur 7; Abidjan TV National-International
## 42788    2.9000             1      Daily Post (Nigeria)               National
##                                                                                                                                                                                                                                                                notes
## 42787 On 31 December 2020, unidentified individuals opened fire on the road toll station of Thomasset town, in the Eastern Ivory Coast, assumedly around Abengourou town (Indenie-Djuablin, Comoe). They destroyed properties. No fatality and/or casualty reported.
## 42788                                                      On 31 December 2020, youth demonstrators set up bonfires and blocked the International Road in Oke Odan (Egbado South LGA, Ogun) over the killing of another youth by Nigeria Customs Service operatives.
##       fatalities                 tags  timestamp
## 42787          0                      1610409617
## 42788          0 crowd size=no report 1610409633

Visualisation de la carte interactive pour l’année 2020

Create_map(AOI_Raster_2020)

## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Tableau du nombre d’attaques par année et par type d’événement

t <- data %>% 
  group_by(year, event_type) %>%
  summarise(attacks_number =n())

## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.

Tableau du nombre total d’attaques par année

t1 <- data %>% 
  group_by(year) %>%
  summarise(attacks_number =n())

Visualisation du nombre total d’attaques par année

ggplot(t1, aes(y= attacks_number, x=year))+
  geom_path(linewidth= 1.2, color="blue")+
  geom_point(size=2, color="red")+
  theme_minimal() +
  labs(title = "Nombre d'attaques par année") 

Visualisation du nombre d’attaques par type d’événement pour chaque année

ggplot(t, aes(y= attacks_number, x=year, color = event_type, group=event_type))+
  geom_path(linewidth= 1.2)+
  geom_point(size=2)+
  theme_minimal() +
  labs(title = "Nombre d'attaques par année")

TP 11 EN LUI MÊME

1- Binarisons le raster evenements

AOI_Raster_binaire <- AOI_Raster_2020
values(AOI_Raster_binaire) <- ifelse(values(AOI_Raster_2020) >= 5, 1, 0) # Les valeurs 1 sont prises si le nombre d'evenements est supérieur à 5

Sauvegarde du raster binaire

writeRaster(AOI_Raster_binaire, "Rasterisation_2020_Binaire_EVENEMENTS.tif", format = "GTiff", overwrite = TRUE)

Verification qu’il n’y a que des valeurs binaires

summary(values(AOI_Raster_binaire))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    0.00    0.00    0.25    0.50    1.00   94959

unique(values(AOI_Raster_binaire))

## [1] NA  1  0

Visualisation

## Palette de couleurs
pal <- colorFactor(palette = c("transparent", "blue"), domain = c(0, 1)) # Couleur bleue si 5evenements et plus
## Visualisation interactive avec Leaflet : Créer une carte interactive pour visualiser le raster binaire
leaflet() %>%
  addTiles() %>%  # Couche de base (OpenStreetMap)
  # Ajouter les limites administratives
  addPolygons(data = shp, color = "brown", weight = 2, opacity = 1, fillOpacity = 0.2,
              popup = ~ADM1_FR) %>%  # Popup avec le nom des régions
  # Ajouter l'image raster binaire
  addRasterImage(AOI_Raster_binaire, colors = pal, opacity = 0.8) %>%
  # Ajouter une légende
  addLegend(pal = pal, values = c(0, 1),
            title = "Valeurs binaires",
            labels = c("Moins de 5 événements", "5 événements ou plus"),
            position = "bottomright") %>%
# Ajouter des contrôles supplémentaires
  addResetMapButton() %>%  # Bouton pour recentrer la carte
  addFullscreenControl()   # Contrôle pour basculer en plein écran

## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

VERIFICATIONS

a. Vérifier les valeurs uniques

valeurs_uniques <- unique(values(AOI_Raster_binaire))
print("Valeurs uniques du raster binaire :" )

## [1] "Valeurs uniques du raster binaire :"

print(valeurs_uniques)

## [1] NA  1  0

b. Faire une table de comparaison pour voir si la binarisation est bien faite

# Créer et binariser directement en excluant les NA
AOI_Raster_binaire[!is.na(values(AOI_Raster)) & values(AOI_Raster) >= 5] <- 1
AOI_Raster_binaire[!is.na(values(AOI_Raster)) & values(AOI_Raster) < 5] <- 0

# Créer la table de comparaison sans NA
comparaison <- data.frame(
  Original = values(AOI_Raster)[!is.na(values(AOI_Raster))],
  Binaire = values(AOI_Raster_binaire)[!is.na(values(AOI_Raster))]
)
# Afficher un échantillon
head(comparaison, 20)

##    Original Binaire
## 1         6       1
## 2         1       0
## 3         2       0
## 4         7       1
## 5         4       0
## 6         2       0
## 7         1       0
## 8         1       0
## 9         3       0
## 10       36       1
## 11       99       1
## 12        6       1
## 13        1       0
## 14       12       1
## 15        2       0
## 16        2       0
## 17       44       1
## 18        3       0
## 19        1       0
## 20        1       0

D’après cet échantillon, on voit que ca a bien fait la binarisation

2- Binarisons le raster population: plus de 50 hbts à 1 sinon 0

Importer le raster de population

pop <- raster("mli_ppp_2020_constrained.tif")

Calcul du facteur d’agrégation (facteur = résolution cible / résolution actuelle)

fact <- round(5000/100)

Agréger le raster à 5 km en utilisant la somme

pop_5km <- aggregate(pop, fact = fact, fun = sum, na.rm = TRUE)

Sauvegarder le raster agrégé

writeRaster(pop_5km, "mli_ppp_2020_constrained_5km.tif", format = "GTiff", overwrite = TRUE)

VERIFICATION

a. Afficher et comparer avant/après

cat("Somme totale avant agrégation :", cellStats(pop, sum, na.rm = TRUE), "\n")

## Somme totale avant agrégation : 23268169

cat("Somme totale après agrégation :", cellStats(pop_5km, sum, na.rm = TRUE), "\n")

## Somme totale après agrégation : 23268169

On remarque que c’est la même chose.

Vérifions que les systèmes de coordonnées sont compatibles

crs(pop_5km) # En degres (CRS=WGS84)

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

crs(AOI_Raster_binaire) # En metre (UTM zone 29N)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 29N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16029]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

Les systèmes de coordonnées ne sont pas les mêmes, on va donc reprojeter

Reprojection avec méthode nearest neigbours

pop_5km_utm <- projectRaster(pop_5km, crs = crs(AOI_Raster_binaire), method = "ngb")
crs(pop_5km_utm) # C'est maintenant en en metre (UTM zone 29N)

## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs 
## WKT2 2019 representation:
## PROJCRS["unknown",
##     BASEGEOGCRS["unknown",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8901]]],
##     CONVERSION["UTM zone 29N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]],
##         ID["EPSG",16029]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

Verification de la correspondance de resolution

res(pop_5km_utm)

## [1] 4440 4620

res(AOI_Raster_binaire)

## [1] 5000 5000

Resample pour qu’on ait également la meme dimension

pop_resampled <- resample(pop_5km_utm, AOI_Raster_binaire, method = "ngb")

Rebinarisation stricte

pop_resampled_binary <- calc(pop_resampled, fun = function(x) {
  ifelse(x < 50, 0, 1)
})

Vérification finale

print(unique(values(pop_resampled_binary)))

## [1] NA  1  0

3- RASTER PRODUIT

Vérifier les dimensions des deux rasters avant de pouvoir faire la multiplication

dim_AOI <- dim(AOI_Raster_binaire)
dim_pop <- dim(pop_resampled_binary)
print(paste("Dimensions du raster AOI_Raster_binaire : ", paste(dim_AOI, collapse = " x ")))

## [1] "Dimensions du raster AOI_Raster_binaire :  277 x 350 x 1"

print(paste("Dimensions du raster pop_5km_binary : ", paste(dim_pop, collapse = " x ")))

## [1] "Dimensions du raster pop_5km_binary :  277 x 350 x 1"

print(unique(values(pop_resampled_binary))) 

## [1] NA  1  0

print(unique(values(AOI_Raster_binaire)))

## [1] NA  1  0

FAISONS ACTUELLEMENT LA MULTIPLICATION

mult_raster <- AOI_Raster_binaire * pop_resampled_binary

Vérifier un résumé des valeurs du raster résultant

summary(values(mult_raster))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    0.00    0.00    0.23    0.00    1.00   95104

4- Calculer par admin, du CDI

AU NIVEAU DES REGIONS

shp_region <- st_read("DONNEES_MALI/mli_admbnda_adm1_1m_gov_20211220.shp", quiet= TRUE)

Calculer le nombre de 1 dans le raster binarisé pour chaque zone administrative

pop_count <- extract(pop_resampled_binary, shp_region, fun = sum, na.rm = TRUE)

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

print(pop_count)

##       [,1]
##  [1,] 3012
##  [2,] 2491
##  [3,] 2590
##  [4,] 1888
##  [5,] 2253
##  [6,] 1384
##  [7,] 1371
##  [8,]  424
##  [9,]    9
## [10,]  274

Calculer le nombre de 1 dans le raster produit pour chaque zone administrative

prod_count <- extract(mult_raster, shp_region, fun = sum, na.rm = TRUE)

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

print(prod_count)

##       [,1]
##  [1,]   14
##  [2,]   14
##  [3,]   12
##  [4,]   73
##  [5,]  192
##  [6,]   30
##  [7,]   70
##  [8,]   10
##  [9,]    1
## [10,]   11

Calculer le CDI pour chaque zone administrative

CDI <- prod_count / pop_count

Afficher les résultats

data.frame(Admin = shp_region$ADM1_FR, CDI = CDI)

##         Admin         CDI
## 1       Kayes 0.004648074
## 2   Koulikoro 0.005620233
## 3     Sikasso 0.004633205
## 4       Ségou 0.038665254
## 5       Mopti 0.085219707
## 6  Tombouctou 0.021676301
## 7         Gao 0.051057622
## 8       Kidal 0.023584906
## 9      Bamako 0.111111111
## 10     Menaka 0.040145985

Visualisation

# Joindre les valeurs du CDI à la carte des régions
shp_region$CDI <- CDI
shp_region_utm <- st_transform(shp_region, crs = 32629)
ggplot(data = shp_region_utm) +
  geom_sf(aes(fill = CDI), color = "white", size = 0.2) +  # Remplir selon CDI
  scale_fill_viridis_c(option = "C", na.value = "gray") +  # Couleur pour les CDI, gérer les NA en mettant le gris
  labs(title = "Carte du CDI par Région", fill = "CDI") +
    geom_sf_text(aes(label = ADM1_FR), size = 3, color = "black", fontface = "bold") +
    annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm")) +
    theme_minimal() +
    theme(legend.position = "right",  
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 8)) +
  guides(fill = guide_colorbar(title = "CDI", title.position = "top", title.hjust = 0.5, barheight = 10))

AU NIVEAU DU PAYS

shp_pays <- st_read("DONNEES_MALI/mli_admbnda_adm0_1m_gov_20211220.shp", quiet= TRUE)
pop_count <- extract(pop_resampled_binary, shp_pays, fun = sum, na.rm = TRUE)

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

print(pop_count)

##       [,1]
## [1,] 15691

prod_count <- extract(mult_raster, shp_pays, fun = sum, na.rm = TRUE)

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

print(prod_count)

##      [,1]
## [1,]  427

CDI <- prod_count / pop_count
data.frame(Admin = shp_pays$ADM0_FR, CDI = CDI)

##       Admin        CDI
## 1 Mali (le) 0.02721305

AU NIVEAU DU DEPARTEMENT

shp_departement <- st_read("DONNEES_MALI/mli_admbnda_adm2_1m_gov_20211220.shp", quiet= TRUE)
pop_count <- extract(pop_resampled_binary, shp_departement, fun = sum, na.rm = TRUE)

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

print(pop_count)

##       [,1]
##  [1,]  481
##  [2,]  366
##  [3,]  470
##  [4,]  400
##  [5,]  868
##  [6,]  293
##  [7,]  134
##  [8,]  221
##  [9,]  468
## [10,]  150
## [11,]  575
## [12,]  337
## [13,]  205
## [14,]  535
## [15,]  696
## [16,]  200
## [17,]  277
## [18,]  381
## [19,]  591
## [20,]  258
## [21,]  187
## [22,]  159
## [23,]  235
## [24,]  228
## [25,]  323
## [26,]  231
## [27,]  466
## [28,]  246
## [29,]  287
## [30,]  236
## [31,]  147
## [32,]  561
## [33,]  336
## [34,]  230
## [35,]  276
## [36,]  180
## [37,]  100
## [38,]  397
## [39,]  464
## [40,]  236
## [41,]  187
## [42,]  528
## [43,]  247
## [44,]  596
## [45,]   74
## [46,]  182
## [47,]   96
## [48,]   72
## [49,]    9
## [50,]   88
## [51,]   90
## [52,]   47
## [53,]   49

prod_count <- extract(mult_raster, shp_departement, fun = sum, na.rm = TRUE)

## Warning in .local(x, y, ...): Transforming SpatialPolygons to the crs of the
## Raster

print(prod_count)

##       [,1]
##  [1,]    1
##  [2,]    4
##  [3,]    3
##  [4,]    1
##  [5,]    2
##  [6,]    2
##  [7,]    1
##  [8,]    3
##  [9,]    1
## [10,]    0
## [11,]    1
## [12,]    2
## [13,]    1
## [14,]    6
## [15,]    1
## [16,]    0
## [17,]    1
## [18,]    1
## [19,]    1
## [20,]    1
## [21,]    7
## [22,]    0
## [23,]    2
## [24,]   23
## [25,]   28
## [26,]    3
## [27,]    8
## [28,]    9
## [29,]   44
## [30,]   27
## [31,]   20
## [32,]   29
## [33,]   34
## [34,]   25
## [35,]   10
## [36,]    3
## [37,]    1
## [38,]    6
## [39,]   10
## [40,]    8
## [41,]    5
## [42,]   37
## [43,]    6
## [44,]   27
## [45,]    2
## [46,]    5
## [47,]    3
## [48,]    0
## [49,]    1
## [50,]    8
## [51,]    2
## [52,]    1
## [53,]    0

CDI <- prod_count / pop_count
data.frame(Admin = shp_departement$ADM2_FR, CDI=CDI)

##             Admin         CDI
## 1       Bafoulabe 0.002079002
## 2           Diéma 0.010928962
## 3           Kayes 0.006382979
## 4         Kéniéba 0.002500000
## 5            Kita 0.002304147
## 6           Nioro 0.006825939
## 7        Yélimané 0.007462687
## 8         Banamba 0.013574661
## 9          Dioila 0.002136752
## 10        Kangaba 0.000000000
## 11           Kati 0.001739130
## 12       Kolokani 0.005934718
## 13      Koulikoro 0.004878049
## 14           Nara 0.011214953
## 15       Bougouni 0.001436782
## 16        Kadiolo 0.000000000
## 17     Kolondieba 0.003610108
## 18       Koutiala 0.002624672
## 19        Sikasso 0.001692047
## 20      Yanfolila 0.003875969
## 21        Yorosso 0.037433155
## 22      Baraouéli 0.000000000
## 23            Bla 0.008510638
## 24         Macina 0.100877193
## 25          Niono 0.086687307
## 26            San 0.012987013
## 27          Ségou 0.017167382
## 28       Tominian 0.036585366
## 29     Bandiagara 0.153310105
## 30        Bankass 0.114406780
## 31         Djenné 0.136054422
## 32       Douentza 0.051693405
## 33           Koro 0.101190476
## 34          Mopti 0.108695652
## 35       Ténenkou 0.036231884
## 36       Youwarou 0.016666667
## 37           Diré 0.010000000
## 38        Goundam 0.015113350
## 39 Gourma-Rharous 0.021551724
## 40       Niafunké 0.033898305
## 41     Tombouctou 0.026737968
## 42        Ansongo 0.070075758
## 43         Bourem 0.024291498
## 44            Gao 0.045302013
## 45       Abeibara 0.027027027
## 46          Kidal 0.027472527
## 47       Tessalit 0.031250000
## 48     Tin-Essako 0.000000000
## 49         Bamako 0.111111111
## 50         Ménaka 0.090909091
## 51 Anderamboukane 0.022222222
## 52         Inekar 0.021276596
## 53      Tidermene 0.000000000